home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / VideoToolbox 96.06.15 / (Utilities) / CalibrateLuminance / CalibrateLuminance.c < prev    next >
Text File  |  1996-03-08  |  44KB  |  1,067 lines
  1. /*
  2. CalibrateLuminance.c
  3. by Denis Pelli, Lan Zhang, and Preeti Verghese
  4. You will need "Numerical Recipes in C" in order to compile this file; see
  5. limitations below.
  6.  
  7. USE
  8.  
  9. You must run this program (or your own equivalent) to calibrate your video card,
  10. ISR Video Attenuator, and monochrome monitor, as described by Pelli & Zhang
  11. (1991).
  12.  
  13. D.G. Pelli and L. Zhang (1991) Accurate control of contrast on microcomputer
  14. displays. Vision Research, 31:1337-1360.
  15.  
  16. The program's results are stored in a file called LuminanceRecord?.h that you
  17. can then use in the rest of your C programs. "?" will be a digit indicating
  18. which monitor. LuminanceRecord?.h describes the gains of the three video
  19. pathways of your video card and ISR Video Attenuator, and describes the gamma
  20. function of your monitor, to allow automatic gamma correction later, using
  21. SetLuminance(), etc.
  22.  
  23. There are two ways of using the LuminanceRecord file. You can use the
  24. preprocessor #include statement, in the midst of your C code, provided you've
  25. got a LuminanceRecord structure called LR. Or you can use the VideoToolbox
  26. routine called ReadLuminanceRecord() to read the data (at runtime) into a
  27. LuminanceRecord of any name.
  28.  
  29. For historical reasons this program not only measures the gamma function, but
  30. also fits a variety of functions to it. In fact the latest implementation of
  31. Luminance.c, which is the only routine that makes direct use of the gamma
  32. description, would be quite content with just the tabulated gamma function. At
  33. present the gamma function is tabulated very coarsely by CalibrateLuminance, so
  34. it doesn't bother to put this table into the header file LuminanceRecord.h.
  35. However, if you'd rather write your own substitute for this program, then I
  36. would suggest just measuring the luminance at 256 different levels (from 0 to
  37. VMax) and saving the table in the header file. You will also need to measure
  38. the RGB gains, but that can easily be done directly, using an oscilloscope,
  39. rather than trying to be clever, as in here, where we invert the nonlinearity to
  40. infer voltages from the nonlinear luminances that we measure.
  41.  
  42. CalibrateLuminance also produces a CricketGraph file, CalibrateLuminance%d.data
  43. (where %d is replaced by the monitor number), which is suitable for graphing the
  44. monitor's gamma function. Use the CricketGraph format file CalibrateLuminance.format.
  45.  
  46. This program measures the luminance of a patch on the screen, using a photometer
  47. and, optionally, a 12-bit Analog to Digital Converter (ADC). If you have a Data
  48. Translation Forerunner Analog-to-Digital card, then GetVoltage.c will
  49. automatically find and use it. You can set the MANUAL flag to force manual
  50. operation. If your photometer has an analog output, it's convenient and more
  51. accurate to have all the voltages read in automatically via an ADC, but you
  52. won't need to run CalibrateLuminance very often, and you can get by with manual
  53. calibration.
  54.  
  55. You should set the background luminance to approximately the same value as you
  56. will use in your experiment. Note that the channels gains depend on the DACs as
  57. well as the ISR Video Attenuator, and will vary from DAC to DAC by ±5%. That's
  58. why you must do this calibration yourself using your own video card and monitor.
  59.  
  60. This program can calibrate any of your screens, including the main screen. It
  61. now works fine with a single screen, alternating the user dialog with the
  62. measurements.
  63.  
  64. LIMITATIONS
  65.  
  66. CalibrateLuminance uses routines from Numerical Recipes in C to do the
  67. polynomial and power law fits to the gamma function. They're copyrighted, so I
  68. can't distribute them. Note: I HAVE included a compiled application
  69. CalibrateLuminance that you can use NOW. You only need to buy the Numerical
  70. Recipes if you want to MODIFY CalibrateLuminance.c. See "Improve Numerical Recipes" 
  71. in the VideoToolbox Notes folder.
  72.  
  73. HISTORY
  74.  
  75. 4/25/89    Preeti and Denis
  76. 6/18/89    Denis added numerical display of clut index and triple grating for first screen.
  77. 8/4/89    Denis replaced SetEntries call by GDSetEntries, and generally updated everything.
  78. 8/18/89    Lan Measure the whole routine(R+G+B and R G B gain) 40 times to eliminate the
  79.         effect of screen luminance drift.
  80. 9/8/89    Lan & Denis generally tidied it up.
  81. 9/8/89    denis: replaced polynomial fit by powerlaw fit
  82. 9/10/89 denis: got gain measurement to work with sufficient accuracy. 
  83.         This involved many small changes. I now
  84.         wait for a second after any large luminance change to allow the
  85.         photometer to settle. I also measure the channel gains at many different
  86.         settings of the other DACs, in order to average out the effect of DAC
  87.         inaccuracies.
  88. 10/30/89 Lan & Denis: introduced the option of manual readings, and made console smaller.
  89. 11/17/89 Lan & Denis: cleaned up for general release.
  90. 11/30/89 Denis: added comment to LuminanceRecord.h explaining power[] parameters.
  91. 3/29/90    dgp    2.15 Updated to use new GetVoltage that looks for ForeRunner card,
  92.             and new GDOpenWindow() that uses CWindowPtr instead of WindowPtr. Introduced
  93.             conditionals so it compiles without errors under MPW C 3.1. 
  94. 4/21/90    dgp    2.20. Fixed the bug in the fixed-power power law fit. Updated whole file
  95.             to be compatible with latest versions of all subroutines. Corrected printout
  96.             of nBackground.
  97. 7/28/90    dgp    2.3. Added code (in GetALuminance.c) to find the equivalent number to
  98.             produce any desired background luminance.
  99.             All luminance measurement is now done via the new subroutine GetALuminance. 
  100.             Automatic and manual measurements now use the exactly the same code.
  101.             Calibration of (linear) RGB gains is now optional.
  102. 9/18/90    dgp    2.4. Changed all instances of "v" to "V". Pelli
  103.             & Zhang (1991) refer to a nominal voltage v; this file
  104.             now refers to the "equivalent number" V; they are related by V=VMax*v.
  105. 9/22/90    dgp    2.5. Fixed cosmetic errors that prevented use of this program on the main
  106.             screen. The trick is appropriate use of BringToFront() and SendBehind()
  107.             referring only to my own window. Trying to do BringToFront() on the console
  108.             often caused a bus error, for no obvious reason. I'll have to ask Mike
  109.             Kahl.
  110. 9/24/90    dgp    2.6. Added screen to the LuminanceRecord.h file and appended it to
  111.             the LuminanceRecord.h and .data filenames. This makes it easy to calibrate
  112.             all your monitors and keep the records straight.
  113.             Added LR.date string.
  114.             New default is to retain old gainAccuracy when rgb gains are not remeasured.
  115. 10/10/90 dgp Added SetDepth(). Fixed bug that reducing the number of frames sampled
  116.             per a/d measurement by the number of cycles. The number of frames per
  117.             measurement is now fixed.
  118. 10/12/90 dgp Added LR.dpi and LR.Hz to the LuminanceRecord?.h file.
  119. 10/17/90 dgp Removed unused variables. Added reference to paper to LuminanceRecord?.h.
  120. 12/12/90 dgp Measure and subtract off the dark voltage.
  121.             Fixed dumb error of assuming wrong type when printing LP->VMax
  122.             & LP->coefficients, which was resulting in zeroes in the LuminanceRecord.
  123. 4/15/91    dgp Check for NewPaletteManager().
  124. 8/24/91    dgp    made compatible with THINK C 5. This required introducing several
  125.             casts: (unsigned long *) and (unsigned char *). Hopefully, this won't
  126.             compromise compatibility with old THINK C 4.
  127. 6/23/92    dgp    added quick check of photometer gain.
  128. 8/27/92    dgp    replace SysEnvirons() by Gestalt()
  129. 10/23/92 dgp read latest LuminanceRecord for default values. Mentioned ReadLuminanceRecord()
  130.             in the LuminanceRecord header file.
  131. 10/24/92 dgp use GDFrameRate() to measure the frame rate.
  132. 11/13/92 dgp advertise ReadLuminanceRecord.c
  133. 12/17/92 dgp Enhanced to support arbitrary dacSize. Get dacSize by calling GDGetGamma(). 
  134.             Mostly I just replaced 255 by LP->VMax. Didn't test it after changes.
  135. 12/21/92 dgp No longer load unused dac bits.
  136. 2/23/93    dgp    Use new GDOpenWindow1 and GDDisposeWindow1.
  137. 4/3/93    dgp    Fix bug that causes GetALuminance to fail (i.e. always set entry to black)
  138.             if the LuminanceRecord hasn't been initialized.
  139. 4/18/93    dgp    Replaced ctSize by clutSize. Now call GDClutSize. Fixed call to SetDepth.
  140. 5/4/93    dgp    Leave gray/color mode alone, but warn if in gray mode.
  141. 5/12/93    dgp When in gray mode don't bother to offer to calibrate rgb gains.
  142. 6/4/93    dgp Conditionally changed the data-saving subroutine to use the
  143.             new WriteLuminanceRecord subroutine. This hasn't been tested yet, and
  144.             we can revert back to the old code by changing the #if 1 to #if 0.
  145.             The virtue of the new code is that it verifies what was written.
  146. 5/28/94 dgp Changed the second argument of FillRect() to be "&qd.black" instead of 
  147. "qd.black" for compatibility with Apple's redefinition of Pattern as a struct 
  148. instead of an array. This makes the call compatible with Apple's Universal Headers.
  149. 5/31/94    dgp    added (ConstPatternParam) cast to above, to retain compatibility
  150. with the old definition of Pattern.
  151. 9/5/94 dgp removed assumption in printf's that int==short.
  152. 11/17/94 dgp updated to use ChooseScreen. Added conditional code for CodeWarrior.
  153. 7/26/95 dgp Josh Solomon reported that it wasn't working. Fixed, mainly by removing the conditional
  154.             CodeWarrior code that was needed for that version, but which is superfluous, and fails,
  155.             under the current version. 
  156.             Added code to restore screen to its original gamma, depth, and mode.
  157. 3/8/96 dgp default of no repeat when !automatic. GDOpenWindow1()'s call to GDGrayColorTable() causes the 
  158.             window frame to become orange. I don't understand why, but it's only a cosmetic bug for this 
  159.             program, so i'm ignoring it for now.
  160. */
  161. #include "VideoToolbox.h"
  162. #include "Luminance.h"
  163. #include "nr.h"            /* prototypes of Numerical Recipes and definition of FLOAT */
  164. #include "nrutil.h"
  165. #if UNIVERSAL_HEADERS
  166.     #include <LowMem.h>
  167. #else
  168.     #define LMGetMBarHeight() (* (short *) 0x0BAA)
  169.     #define LMSetMBarHeight(mBarHeight) ((* (short *) 0x0BAA) = (mBarHeight))
  170. #endif
  171. long InitializeLuminanceRecord(LuminanceRecord *LP,short flags);
  172.  
  173. #define AUTOMATIC        1    /*     =0: always request manual photometer readings
  174.                                    =1: if possible, use A/D converter to read from the photometer directly */
  175. #define NLRGB            16    /* number of times to measure each RGB gain, should be >1 */
  176. #define LUMINANCES        32    /* approx. number of different luminances to measure */
  177. #define    NFRAME            60    /* number of frames per ADC luminance measurement */
  178. #define ROUND_LUMINANCES (2+(256+256/LUMINANCES-1)/(256/LUMINANCES))
  179.                             /* round LUMINANCES up to the next divisor of 256, add 2 */
  180.  
  181. void main(void);
  182. void CalibrateLuminance(void);
  183. void SaveData(LuminanceRecord *LP,int nL,int V[],double L[]
  184.     ,int nLrgb,int nrgb[3][NLRGB][3],double Lrgb[3][NLRGB],double gain[3]);
  185. void SaveLuminanceRecord(LuminanceRecord *LP);
  186. FLOAT PowerRMSError(FLOAT p[]);
  187. double GetALuminance(LuminanceRecord *LP,GDHandle device,
  188.     int frames,double LuminancePerVoltage,int entry,int red,int green,int blue);
  189. void LToVHunt(LuminanceRecord *LP,GDHandle device,CWindowPtr window,
  190.     double LuminancePerVoltage,int frames,double darkLuminance);
  191.  
  192. /* These variables are out here, as globals, so that PowerRMSError can access them directly */
  193. static int nL=ROUND_LUMINANCES;
  194. static double L[ROUND_LUMINANCES];    /* luminance at V[i] */
  195. static int V[ROUND_LUMINANCES];
  196. static int variables=4;
  197. static FLOAT *p;
  198.  
  199. void CalibrateLuminance(void);
  200.  
  201. void main(void)
  202. {
  203.     long system;
  204.     
  205.     Gestalt(gestaltSystemVersion,&system);
  206.     if(system<0x605)PrintfExit("Sorry, I need at least System 6.05.\n");
  207.     Require(gestalt8BitQD);
  208.     StackGrow(sizeof(LuminanceRecord)+2L*3*NLRGB*sizeof(double)+1000);
  209.     CalibrateLuminance();
  210. }
  211.  
  212. void CalibrateLuminance(void)
  213. {
  214.     int i,j,k,ii,bottom,screen,error;
  215.     short fontNumber;
  216.     char string[1000];
  217.     Rect rect,testRect,labelRect;
  218.     GDHandle device=NULL,oldDevice=NULL;
  219.     WindowPtr window=NULL,oldWindow=NULL;
  220.     int nLrgb=NLRGB;
  221.     double Lrgb[3][NLRGB];
  222.     int nrgb[3][NLRGB][3];
  223.     double gain[3],luminancePerVoltage=1000.0;
  224.     LuminanceRecord LR,*LP;
  225.     FLOAT *x,*y,*sig;                        // for polynomial fit
  226.     FLOAT *a,**u,**v,*w,**cvm,chisq;        // for polynomial fit
  227.     int ma;                                    // for polynomial fit
  228.     FLOAT **xi,ftol,fret;                    // for power law fit
  229.     int iter;                                // for power law fit
  230.     double e,f,VV,luminance,darkLuminance;
  231.     int readNumber,readTotal,isGray,testSize,cycles,frames,clutSize,mode;
  232.     Boolean automatic,skipRGB,yes;
  233.     char colorGray[2][8]={"color","gray"};
  234.     short oldPixelSize,oldIsColor;
  235.     
  236.     assert(StackSpace()>4000);
  237.     automatic=AUTOMATIC && (CardSlot(".ForeRunner")>0);    /* disable if card is missing */
  238.     LP=&LR;
  239.     #define WIDTH 105
  240.     #if (THINK_C || THINK_CPLUS || SYMANTEC_C)
  241.         console_options.title="\pCalibrateLuminance";    /* change the console window */
  242.         console_options.nrows=12;
  243.         console_options.top=20;
  244.         console_options.left=1;
  245.         console_options.ncols=WIDTH;
  246.         console_options.txSize=9;
  247.         printf("\n");                                    /* ask THINK C to initialize QuickDraw */
  248.         CopyQuickDrawGlobals();                            // Make sure qd is valid.
  249.     #elif __MWERKS__
  250.         SIOUXSettings.toppixel=LMGetMBarHeight()+1;    // allow for menu bar only
  251.         SIOUXSettings.leftpixel=1;
  252.         SIOUXSettings.rows=13;
  253.         SIOUXSettings.columns=WIDTH;
  254.         SIOUXSettings.autocloseonquit=1;
  255.         SIOUXSettings.showstatusline=0;
  256.         SIOUXSettings.asktosaveonclose=0;
  257.         SIOUXSettings.fontsize=9;
  258.         printf("\n");
  259.         SIOUXSetTitle("\pCalibrateLuminance");
  260.     #else
  261.         InitGraf(&qd.thePort);
  262.         InitWindows();
  263.         InitFonts();
  264.     #endif
  265.  
  266.     printf("\n");    // make sure that oldWindow is the console
  267.     GetPort(&oldWindow);
  268.     oldDevice=GetGDevice();
  269.     SetPort(oldWindow);
  270.     SetGDevice(oldDevice);
  271.  
  272.     printf(BreakLines("Welcome to CalibrateLuminance version %s.\n"
  273.     "This program helps you to calibrate your ISR Video Attenuator, as well as your video card "
  274.     "and monitor. 'Enter' indicates that you should type in a number or letter followed by 'RETURN'. "
  275.     "Most questions provide a default answer (in parentheses). "
  276.     "Just hit 'RETURN' if you want the default. (Or hit Command-. to quit.) "
  277.     "Where appropriate, the default is taken from latest LuminanceRecord.\n",WIDTH),__DATE__);
  278.  
  279.     if(GetScreenDevice(1) == NULL)screen=0;
  280.     else do{
  281.         screen=1;
  282.         screen=ChooseScreen(screen,"Which screen do you want to calibrate?");
  283.     }while(GetScreenDevice(screen)==NULL);
  284.     sprintf(string,"LuminanceRecord%d.h",screen);
  285.     InitializeLuminanceRecord(LP,0);
  286.     i=ReadLuminanceRecord(string,LP,0);        // try to read latest LuminanceRecord
  287.     if(i>0)printf("Successfully read old “%s” to set default values.\n",string);
  288.     else{
  289.         printf("Couldn't find an old copy of “%s”; creating new one from scratch.\n",string);
  290.         LR.dpi=76.0;    /* pixels per inch */
  291.         LP->units="cd/m^2";
  292.     }
  293.     LP->screen=screen;
  294.     LP->rangeSet=0;    /* indicate that range parameters have yet to be set */
  295.     LP->L.exists=0;    /* indicate that luminance table has yet to be initialized */
  296.     GetPort(&oldWindow);
  297.     oldDevice=GetGDevice();
  298.  
  299.     device=GetScreenDevice(LP->screen);        // choose which screen to test
  300.     oldPixelSize=(**(**device).gdPMap).pixelSize;
  301.     oldIsColor=TestDeviceAttribute(device,gdDevType);
  302.     GDSaveGamma(device);
  303.     window=GDOpenWindow1(device);
  304.     SendBehind(window,NULL);
  305.     LP->dacSize=GDDacSize(device);            // Takes 200 µs.
  306.     clutSize=GDClutSize(device);    
  307.     LP->VMin=0;
  308.     LP->VMax=(1L<<LP->dacSize)-1;    /* maximum value that can be loaded into DAC */
  309.     isGray=!TestDeviceAttribute(device,gdDevType);
  310.     printf("Test screen is set to %d colors and %s mode.\n",clutSize,colorGray[isGray]);
  311.     if(clutSize<4){
  312.         for(i=4;i<=32;i*=2){
  313.             mode=HasDepth(device,i,0,0);
  314.             if(mode!=0)break;
  315.         }
  316.         if(mode==0)PrintfExit("\n\007Sorry. This program cannot calibrate a display"
  317.             "with less than 4 colors/grays.");
  318.         printf("Changing screen depth to %d-bit pixels\n",i);
  319.         if(mode==0x100)PrintfExit(BreakLines("Sorry, the HasDepth routine doesn't work properly. "
  320.         "Use the Monitors Control Panel to set the colors/grays and then run "
  321.         "CalibrateLuminance again.\n",WIDTH));
  322.         SetDepth(device,mode,0,0);
  323.         clutSize=GDClutSize(device);    
  324.     }
  325.     if(isGray){
  326.         yes=Choose(0,"Would you prefer color mode?\n",noYes,2);
  327.         if(yes){
  328.             if(HasDepth(device,(**(**device).gdPMap).pixelSize,1,1))
  329.                 SetDepth(device,(**device).gdMode,1,1);
  330.             else printf("Sorry, there is no color mode. Continuing in gray mode.\n");
  331.         }
  332.         isGray=!TestDeviceAttribute(device,gdDevType);
  333.     }else{
  334.         printf("Hit return when ready to continue.\n");
  335.         gets(string);
  336.     }
  337.     LP->Hz=GDFrameRate(device);
  338.     
  339.     if(isGray)skipRGB=1;
  340.     else{
  341.         printf(BreakLines("\nThere are two kinds of measurement:\n"
  342.         "1. Compulsory calibration of the (nonlinear) monitor's gamma function. "
  343.         "This will be affected by your choice of background luminance.\n"
  344.         "2. Optional calibration of the (linear) ISR Video Attenuator 3 channel gains. "
  345.         "This is independent of the background. It takes about 15 minutes.\n\n",WIDTH));
  346.         skipRGB=Choose(0,"Do you wish to skip 2., the optional three-channel calibration?\n",noYes,2);
  347.     }
  348.     if(skipRGB)printf("Skipping 3 channel calibration.\n");
  349.     else printf("All 3 channels will be calibrated.\n");
  350.     
  351.     frames=0;
  352.     if(automatic){
  353.         frames=NFRAME;
  354.         printf("Taking photometer readings automatically. "
  355.             "Assuming %.3f milliVolts per %s.\n",1000.0/luminancePerVoltage,LP->units);
  356.     }
  357.     else printf("Please enter the photometer readings,in %s, as prompted.\n",LP->units); 
  358.     printf("\n");
  359.  
  360.     printf("First let's check the photometer's gain setting.\n"
  361.         "Point the photometer at the screen and hit return:");
  362.     gets(string);
  363.     clutSize=GDClutSize(device);
  364.     do{
  365.         printf("\r");
  366.         SetPort(window);
  367.         BringToFront(window);
  368.         FillRect(&window->portRect,(ConstPatternParam)&qd.black);
  369.         darkLuminance=GetALuminance(LP,device,frames,luminancePerVoltage,clutSize-1,0,0,0);
  370.         FillRect(&window->portRect,(ConstPatternParam)&qd.white);
  371.         luminance=GetALuminance(LP,device,frames,luminancePerVoltage,0,LP->VMax,LP->VMax,LP->VMax);
  372.         luminance-=darkLuminance;
  373.         SendBehind(window,NULL);
  374.         printf("\rWhite seems to be %.1f %s brighter than black.\n",luminance,LP->units);
  375.         yes=Choose(0,"Do you want to try that again?",noYes,2);
  376.     }while(yes);
  377.     printf("\n");
  378.  
  379.     printf("Please occlude the photometer so that I can read the dark level.\n"
  380.     "Hit return when ready:");
  381.     gets(string);
  382.     darkLuminance=GetALuminance(LP,device,frames,luminancePerVoltage,0,LP->VMax,LP->VMax,LP->VMax);
  383.     printf("The dark level of %.2f cd/m^2 will be subtracted from all readings\n"
  384.         ,darkLuminance);
  385.     
  386.     testSize=160;
  387.     printf("Enter, in pixels, the diameter of the test patch (%d):",testSize);
  388.     gets(string);
  389.     sscanf((char *)string,"%d",&testSize);
  390.  
  391.     /* draw circle for placement of photocell */
  392.     SetPort(window);
  393.     rect=window->portRect;
  394.     InsetRect(&rect,(rect.right-testSize)/2,(rect.bottom-testSize)/2);
  395.     PmForeColor(255);
  396.     PenSize(10,10);
  397.     FrameOval(&rect);
  398.     SetPort(oldWindow);
  399.     
  400.     printf("Please remove the occluder. \n");
  401.     printf("Place your photometer so as to read the luminance of the "
  402.         "center of the test screen.\n"); 
  403.  
  404.     if(skipRGB || !automatic)cycles=1;
  405.     else cycles=4;
  406.     printf("Then enter number of times you wish to repeat the whole cycle of measurement. (%d):"
  407.         ,cycles);
  408.     gets((char *)string);
  409.     sscanf((char *)string,"%d",&cycles);
  410.     if(cycles<=0)return;
  411.  
  412.     SetPort(window);
  413.     BringToFront(window);
  414.     FillRect(&window->portRect,(ConstPatternParam)&qd.gray);
  415.     luminance=GetALuminance(LP,device,frames,luminancePerVoltage,0,LP->VMax,LP->VMax,LP->VMax);
  416.     luminance-=darkLuminance;
  417.     printf("50%% dithered gray is %.1f %s\n",luminance,LP->units);
  418.     FillRect(&window->portRect,(ConstPatternParam)&qd.white);
  419.     luminance=GetALuminance(LP,device,frames,luminancePerVoltage,0,LP->VMax,LP->VMax,LP->VMax);
  420.     luminance-=darkLuminance;
  421.     printf("White is %.1f %s\n",luminance,LP->units);
  422.     FillRect(&window->portRect,(ConstPatternParam)&qd.black);
  423.     luminance=GetALuminance(LP,device,frames,luminancePerVoltage,0,LP->VMax,LP->VMax,LP->VMax);
  424.     luminance-=darkLuminance;
  425.     printf("Black is %.1f %s\n",luminance,LP->units);
  426.     SendBehind(window,NULL);
  427.  
  428.     sprintf(string,"It is essential that the calibration be made with the same background "
  429.     "luminance as will ultimately be used in your experiments. You may specify the "
  430.     "background luminance either as a luminance in %s, or as a corresponding "
  431.     "DAC value (0..%d).\n",LP->units,(int)LP->VMax);
  432.     printf("%s",BreakLines(string,WIDTH));
  433.     yes=Choose(1,"Do you wish to enter a luminance?\n",noYes,2);
  434.     if(yes){
  435.         printf("Enter the desired background luminance in %s (%f):",LP->units,LP->LBackground);
  436.         gets(string);
  437.         sscanf(string,"%lf",&LP->LBackground);
  438.         printf("Now hunting for the DAC value . . .\n");
  439.         LToVHunt(LP,device,(CWindowPtr)window,luminancePerVoltage,frames,darkLuminance);
  440.         printf("Will use the nearest DAC value, which is %hd\n",LP->VBackground);
  441.     }else{
  442.         printf("Enter the DAC value for the background (%hd):",LP->VBackground);
  443.         gets(string);
  444.         sscanf((char *)string,"%hd",&LP->VBackground);
  445.     }
  446.     LP->table[2].rgb.red=LP->table[2].rgb.green=LP->table[2].rgb.blue=
  447.         LP->VBackground<<LP->leftShift;
  448.     LoadLuminances(device,LP,2,2);
  449.     SetPort(window);
  450.     BringToFront(window);
  451.     PmBackColor(2);
  452.     EraseRect(&window->portRect);
  453.     luminance=GetALuminance(LP,device,frames,luminancePerVoltage,0,LP->VMax,LP->VMax,LP->VMax);
  454.     luminance-=darkLuminance;
  455.     SendBehind(window,NULL);
  456.     SetPort(oldWindow);
  457.     printf("Your background luminance is %.3f %s\n",luminance,LP->units);
  458.  
  459.     /* plan the luminance and gain measurements */
  460.     if(skipRGB)nLrgb=0;
  461.     readTotal=cycles*(ROUND_LUMINANCES+3*nLrgb);
  462.     printf("Now taking %d readings . . .\n",readTotal);
  463.     for (i=0;i<nL;i++) L[i]=0.0;
  464.     k=0;
  465.     for (i=0;i<nL;i++){
  466.         V[i]=k;
  467.         if(k==LP->VMax) {
  468.             if(i+1<nL) V[++i]=k;    /* put two points at VMax, because it's important */
  469.             nL=i+1;
  470.             break;
  471.         }
  472.         k += (LP->VMax+1)/(nL-2);
  473.         if(k>LP->VMax)k=LP->VMax;
  474.     }
  475.     for (k=0;k<3;k++)for(i=0;i<nLrgb;i+=2){
  476.         for(j=0;j<3;j++) nrgb[k][i+1][j]=nrgb[k][i][j]=(LP->VMax+1)*(0.667+0.333*i/nLrgb);
  477.         nrgb[k][i+1][k]=LP->VMax;
  478.         nrgb[k][i][k]=0;
  479.         if(k==2) nrgb[k][i][k]=(LP->VMax+1)/2;    /* take smaller step, because blue gain is highest */
  480.     }
  481.     for(k=0;k<3;k++)for (i=0;i<nLrgb;i++)Lrgb[k][i]=0.0;
  482.     
  483.     /* measure several times to minimize the effect of luminance drift */
  484.     SetPort(window);
  485.     BringToFront(window);
  486.     GetFNum("\pMonaco",&fontNumber);
  487.     TextFont(fontNumber);
  488.     TextSize(36);
  489.     rect=window->portRect;
  490.     bottom=rect.bottom;
  491.     SetRect(&testRect, (rect.right+rect.left-testSize)/2
  492.         ,(rect.top+rect.bottom-testSize)/2
  493.         ,(rect.right+rect.left+testSize)/2
  494.         ,(rect.top+rect.bottom+testSize)/2 );
  495.     
  496.     /* set up background */
  497.     LP->table[2].rgb.red=LP->table[2].rgb.green=LP->table[2].rgb.blue=
  498.         LP->VBackground<<LP->leftShift;
  499.     LoadLuminances(device,LP,1,2);    /* load background into clut */
  500.     PmBackColor(2);                /* our background */
  501.     EraseRect(&rect);
  502.     PmBackColor(1);                /* our test luminance */
  503.     EraseOval(&testRect);
  504.     FlushEvents(everyEvent,0);
  505.     readNumber=0;    /* initialize to count the number of reading */
  506.     for(j=0;j<cycles;j++)    {
  507.         for (k=-1;k<3;k++)    {
  508.             if(k<0) ii=nL;
  509.             else ii=nLrgb;
  510.             for (i=0;i<ii;i++){
  511.                 /* set up test */
  512.                 SetRect(&labelRect,0,bottom-40,300,bottom);
  513.                 SetPort(window);
  514.                 PmBackColor(2);
  515.                 EraseRect(&labelRect);
  516.                 PmForeColor(clutSize-1);                /* black */
  517.                 MoveTo(0,bottom);
  518.                 if(k<0) DrawPrintf("%3d",V[i]);
  519.                 else DrawPrintf("%3d%4d%4d",k,nrgb[k][i][0],nrgb[k][i][1],nrgb[k][i][2]);
  520.                 SetPort(oldWindow);
  521.                 readNumber++;
  522.                 if(!automatic){
  523.                     printf("%d out of %d readings:\t",readNumber,readTotal);    
  524.                 }
  525.                 if(k<0) {
  526.                     luminance=GetALuminance(LP,device,frames,luminancePerVoltage,1
  527.                         ,V[i],V[i],V[i]);
  528.                     luminance-=darkLuminance;
  529.                     L[i] += luminance/cycles;
  530.                 }
  531.                 else {
  532.                     luminance=GetALuminance(LP,device,frames,luminancePerVoltage,1
  533.                         ,nrgb[k][i][0],nrgb[k][i][1],nrgb[k][i][2]);
  534.                     luminance-=darkLuminance;
  535.                     Lrgb[k][i] += luminance/cycles;
  536.                 }
  537.             }
  538.         }
  539.     }
  540.     GDDisposeWindow1(window);
  541.     SetPort(oldWindow);
  542.     SetGDevice(oldDevice);
  543.     GDRestoreGamma(device);
  544.     GDRestoreDeviceClut(device);
  545.     error=SetDepth(device,oldPixelSize,1<<gdDevType,oldIsColor);
  546.     
  547.     /* polynomial and quadratic fits */
  548.     ma=MAX_COEFFICIENTS;
  549.     x=vector(1,nL);
  550.     y=vector(1,nL);
  551.     sig=vector(1,nL);
  552.     a=vector(1,nL);
  553.     u=matrix(1,nL,1,ma);
  554.     v=matrix(1,ma,1,ma);
  555.     w=vector(1,ma);
  556.     cvm=matrix(1,ma,1,ma);
  557.     for(i=0;i<nL;i++){
  558.         x[i+1]=V[i];
  559.         y[i+1]=L[i];
  560.         sig[i+1]=10.0;
  561.     }
  562.     svdfit(x,y,sig,nL,a,ma,u,v,w,&chisq,fpoly);    /* ma-1th order polynomial curve fit */
  563.     svdvar(v,ma,w,cvm);
  564.     if(ma>MAX_COEFFICIENTS)PrintfExit("Error: too many coefficients\007\n");
  565.     LP->coefficients=ma;
  566.     for(i=0;i<ma;i++) LP->p[i]=a[i+1];
  567.     printf("L(V) =");
  568.     for(i=0;i<ma;i++) printf(" + %6g V^%d",a[i+1],i);
  569.     printf(".  chisq %g\n",chisq);
  570.     svdfit(x,y,sig,nL,a,3,u,v,w,&chisq,fpoly);    /* 2nd order polynomial curve fit */
  571.     svdvar(v,3,w,cvm);
  572.     printf("\nquadratic fit:\n\n");
  573.     for(i=0;i<3;i++) LP->q[i]=a[i+1];
  574.     printf("L(V) =");
  575.     for(i=0;i<3;i++) printf(" + %g V^%d",a[i+1],i);
  576.     printf(".  chisq %g\n",chisq);
  577.     free_vector(x,1,nL);
  578.     free_vector(y,1,nL);
  579.     free_vector(sig,1,nL);
  580.     free_vector(a,1,nL);
  581.     free_matrix(u,1,nL,1,ma);
  582.     free_matrix(v,1,ma,1,ma);
  583.     free_vector(w,1,ma);
  584.     free_matrix(cvm,1,ma,1,ma);
  585.     e=0.0;
  586.     for(i=0;i<nL;i++){
  587.         f=0.0;
  588.         VV=1.0;
  589.         for(j=0;j<LP->coefficients;j++){
  590.             f+=LP->p[j]*VV;
  591.             VV*=V[i];
  592.         }
  593.         e+=(L[i]-f)*(L[i]-f);
  594.     }
  595.     LP->polynomialError=sqrt(e/nL);
  596.     e=0.0;
  597.     for(i=0;i<nL;i++){
  598.         f=0.0;
  599.         VV=1.0;
  600.         for(j=0;j<3;j++){
  601.             f+=LP->q[j]*VV;
  602.             VV*=V[i];
  603.         }
  604.         e+=(L[i]-f)*(L[i]-f);
  605.     }
  606.     LP->quadraticError=sqrt(e/nL);
  607.     
  608.     /* power law fit */
  609.     /* L=p[1]+Rectify(p[2]+p[3]*V)^p[4] */
  610.     p=vector(1,4);    /* initial starting point */
  611.     /*
  612.         It is necessary to have a reasonable starting point or the search can get
  613.         stuck out in the boondocks. Since the search is quite slow it is desirable to
  614.         give it as good a starting point as possible. Therefore I make use of the
  615.         quadratic fit to try to give it a pretty good starting point.
  616.     */
  617.     /* use the quadratic fit as the starting point */
  618.     p[1]=LP->q[0]-0.25*LP->q[1]*LP->q[1]/LP->q[2];
  619.     p[2]=0.5*LP->q[1]/sqrt(LP->q[2]);
  620.     p[3]=sqrt(LP->q[2]);
  621.     p[4]=2.0;
  622.     /*
  623.         On second thought, try starting with a gamma of 2.3.
  624.         Set the luminance offset p[1] to the measured dc level.
  625.         Set the brightness p[2] so as to match the number offset in the quadratic fit.
  626.         Set the contrast p[3] so as to fit the maximum data point.
  627.         Set gamma p[4] to 2.3
  628.     */
  629.     p[1]=L[1];
  630.     p[2]=0.5*LP->q[1]/sqrt(LP->q[2]);
  631.     p[4]=2.3;
  632.     p[3]=(pow(L[nL-1]-p[1],1.0/p[4])-p[2])/V[nL-1];
  633.     p[2]=p[3]*0.5*LP->q[1]/LP->q[2];
  634.     p[3]=(pow(L[nL-1]-p[1],1.0/p[4])-p[2])/V[nL-1];
  635.     xi=matrix(1,4,1,4);    /* initial set of directions */
  636.     for(i=1;i<=4;i++)for(j=1;j<=4;j++)xi[i][j]=0.0;
  637.     xi[1][1]=0.1;
  638.     xi[2][2]=1.0;
  639.     xi[3][3]=0.1;
  640.     xi[4][4]=0.01;
  641.     ftol=1e-2;    /* fractional tolerance on function value when done */
  642.     printf("\npower law fit:\n\n");
  643.     variables=4;    /* inform PowerRMSError */
  644.     fret=PowerRMSError(p);
  645.     printf("L(V) = %9.4f + Rectify(%9.4f + %9.4f V)^%9.4f ± %9.4f\n",p[1],p[2],p[3],p[4],fret);
  646.     powell(p,xi,4,ftol,&iter,&fret,&PowerRMSError);
  647.     printf("L(V) = %9.4f + Rectify(%9.4f + %9.4f V)^%9.4f ± %9.4f\n",p[1],p[2],p[3],p[4],fret);
  648.     for(i=0;i<4;i++) LP->power[i]=p[i+1];
  649.     free_vector(p,1,4);
  650.     free_matrix(xi,1,4,1,4);
  651.     e=0.0;
  652.     for(i=0;i<nL;i++){
  653.         f=LP->power[1]+LP->power[2]*V[i];
  654.         if(f>0.0) f=LP->power[0]+pow(f,LP->power[3]);
  655.         else f=LP->power[0];
  656.         e+=(L[i]-f)*(L[i]-f);
  657.     }
  658.     LP->powerError=sqrt(e/nL);
  659.     LP->VMin=0;                    /* minimum value that can be loaded into DAC */
  660.     LP->VMax=(1L<<LP->dacSize)-1;    /* maximum value that can be loaded into DAC */
  661.     LP->LMin=VToL(LP,LP->VMin);    /* min luminance */
  662.     LP->LMax=VToL(LP,LP->VMax);    /* max luminance */
  663.     
  664.     /* power law fit with a FIXED gamma */
  665.     /* This is solely for study of the effects of the contrast and brightness knobs */
  666.     /* L=p[1]+Rectify(p[2]+p[3]*V)^p[4] */
  667.     p=vector(1,4);    /* initial starting point */
  668.     /*
  669.         It is necessary to have a reasonable starting point or the search can get
  670.         stuck out in the boondocks. Since the search is quite slow it is desirable to
  671.         give it as good a starting point as possible. 
  672.     */
  673.     for(i=0;i<4;i++)p[i+1]=LP->power[i];
  674.     p[4]=2.28;
  675.     xi=matrix(1,3,1,3);    /* initial set of directions */
  676.     for(i=1;i<=3;i++)for(j=1;j<=3;j++)xi[i][j]=0.0;
  677.     xi[1][1]=0.1;
  678.     xi[2][2]=1.0;
  679.     xi[3][3]=0.1;
  680.     ftol=1e-2;        /* fractional tolerance on function value when done */
  681.     variables=3;    /* global, to inform PowerRMSError() */
  682.     printf("\npower law fit, with fixed gamma:\n\n");
  683.     fret=PowerRMSError(p);
  684.     printf("L(V) = %9.4f + Rectify(%9.4f + %9.4f V)^%9.4f ± %9.4f\n",p[1],p[2],p[3],p[4],fret);
  685.     powell(p,xi,3,ftol,&iter,&fret,&PowerRMSError);
  686.     printf("L(V) = %9.4f + Rectify(%9.4f + %9.4f V)^%9.4f ± %9.4f\n",p[1],p[2],p[3],p[4],fret);
  687.     for(i=0;i<4;i++) LP->fixedPower[i]=p[i+1];
  688.     free_vector(p,1,4);
  689.     free_matrix(xi,1,3,1,3);
  690.     e=0.0;
  691.     for(i=0;i<nL;i++){
  692.         f=LP->fixedPower[1]+LP->fixedPower[2]*V[i];
  693.         if(f>0.0) f=LP->fixedPower[0]+pow(f,LP->fixedPower[3]);
  694.         else f=LP->fixedPower[0];
  695.         e += (L[i]-f)*(L[i]-f);
  696.     }
  697.     LP->fixedPowerError=sqrt(e/nL);
  698.     
  699.     if(skipRGB){
  700.         if(isGray){
  701.             LP->r=LP->b=0.0;
  702.             LP->g=1.0;
  703.         }else{
  704.             printf("Please enter gain of red channel   (%6.4f):",LP->r);
  705.             gets(string);
  706.             sscanf((char *)string,"%lf",&LP->r);
  707.             printf("Please enter gain of green channel (%6.4f):",LP->g);
  708.             gets(string);
  709.             sscanf((char *)string,"%lf",&LP->g);
  710.             printf("Please enter gain of blue channel  (%6.4f):",LP->b);
  711.             gets(string);
  712.             sscanf((char *)string,"%lf",&LP->b);
  713.         }
  714.         gain[0]=LP->r;
  715.         gain[1]=LP->g;
  716.         gain[2]=LP->b;
  717.     }
  718.     else {    
  719.         /* calculate R,G,B gains */
  720.         for(k=0;k<3;k++){
  721.             gain[k]=0.0;
  722.             for(i=0;i<nLrgb;i+=2)
  723.                 gain[k]+=(LToV(LP,Lrgb[k][i])-LToV(LP,Lrgb[k][i+1]))/(nrgb[k][i][k]-nrgb[k][i+1][k]);
  724.             gain[k] /= nLrgb/2;
  725.             printf("gain[%d]=%6.4f\n",k,gain[k]);
  726.         }
  727.     }
  728.     printf("Sum of gains %6.4f\n",gain[0]+gain[1]+gain[2]);
  729.     LP->r=gain[0]/(gain[0]+gain[1]+gain[2]);        /* Normalize gain */
  730.     LP->g=gain[1]/(gain[0]+gain[1]+gain[2]);
  731.     LP->b=gain[2]/(gain[0]+gain[1]+gain[2]);
  732.     e=1.0-(gain[0]+gain[1]+gain[2]);
  733.     if(skipRGB && fabs(e)<1e-6){
  734.         printf("Please enter gain accuracy   (%6.4f):",LP->gainAccuracy);
  735.         gets(string);
  736.         sscanf((char *)string,"%lf",&LP->gainAccuracy);
  737.     }
  738.     else LP->gainAccuracy=e;
  739.     LP->LBackground=VToL(LP,LP->VBackground);
  740.     LP->gm=(LP->VMax/(2.0*LP->LBackground))
  741.         *(LP->LBackground-VToL(LP,LToV(LP,LP->LBackground)-1.0));
  742.     
  743.     SaveData(LP,nL,V,L,nLrgb,nrgb,Lrgb,gain);
  744.     SaveLuminanceRecord(LP);
  745.     sprintf(string,
  746.         "Congratulations! You've finished the luminance calibration of screen %hd. The calibration results "
  747.         "have been saved as CalibrateLuminance%hd.data and LuminanceRecord%hd.h.\n"
  748.         "     The CalibrateLuminance%hd.data file is a text file suitable for graphing the monitor's gamma function. "
  749.         "We recommend using CricketGraph with the CalibrateLuminance.format. The V column is the equivalent "
  750.         "number loaded into the clut (V=%hd*v), and the L column is the measured luminance in %s.\n"
  751.         "     The program's results have also been stored in a C header file called LuminanceRecord%hd.h that all "
  752.         "your programs should read by using either #include at compile time, or ReadLuminanceRecord() at "
  753.         "runtime. LuminanceRecord%hd.h describes the gains of the three video pathways of your video card and "
  754.         "ISR Video Attenuator, and describes the gamma function of your monitor, to allow automatic gamma "
  755.         "correction later, using SetLuminance(), etc. See Pelli & Zhang (1991) Vision Research 31:1337-1360. Good luck.\n"
  756.         "Hit RETURN to exit:"
  757.         ,LP->screen,LP->screen,LP->screen,LP->screen,LP->VMax,LP->units,LP->screen,LP->screen
  758.     );
  759.     printf("%s",BreakLines(string,WIDTH));
  760.     #if __MWERKS__
  761.         gets(string);
  762.         SIOUXSettings.autocloseonquit=1;
  763.     #endif
  764. }
  765.  
  766. void SaveData(LuminanceRecord *LP,int nL,int V[],double L[]
  767.     ,int nLrgb,int nrgb[3][NLRGB][3],double Lrgb[3][NLRGB],double gain[3])
  768. /* create Cricket data file, to graph the gamma function */
  769. {
  770.     FILE *myfile;
  771.     int i,j,k;
  772.     double VV,LFit;
  773.     static char fileName[80];
  774.     
  775.     sprintf(fileName,"LuminanceCalibration%hd.data",LP->screen);
  776.     myfile=fopen(fileName,"w");
  777.     SetFileInfo(fileName,'TEXT','CGRF');
  778.     fprintf(myfile,"*\n");
  779.     fprintf(myfile,"V\tL \tV(L) ");
  780.     fprintf(myfile,"\tquadratic\tpolynomial\tpower ");
  781.     fprintf(myfile,"\tL-quadratic\tL-polynomial\tL-power ");
  782.     fprintf(myfile,"\tVr \tVg \tVb ");
  783.     fprintf(myfile,"\tLr \tLg \tLb \tV(Lr) \tV(Lg) \tV(Lb) ");
  784.     fprintf(myfile,"\tLR.r\tLR.g \tLR.b \tsum of gains");
  785.     for(i=0;i<nL;i++)    {
  786.         fprintf(myfile,"\n%d\t",V[i]);
  787.         fprintf(myfile,"%9.4g\t",L[i]);
  788.         fprintf(myfile,"%9.4g\t",LToV(LP,L[i]));
  789.         /* quadratic fit */
  790.         LFit=0.0;
  791.         VV=1.0;
  792.         for(j=0;j<3;j++){
  793.             LFit += LP->q[j]*VV;
  794.             VV *= V[i];
  795.         }
  796.         fprintf(myfile,"%9.4g\t",LFit);
  797.         /* polynomial fit */
  798.         LFit=0.0;
  799.         VV=1.0;
  800.         for(j=0;j<LP->coefficients;j++){
  801.             LFit += LP->p[j]*VV;
  802.             VV *= V[i];
  803.         }
  804.         fprintf(myfile,"%9.4g\t",LFit);
  805.         /* power law fit */
  806.         LFit=LP->power[1]+LP->power[2]*V[i];
  807.         if(LFit>0.0) LFit=LP->power[0]+pow(LFit,LP->power[3]);
  808.         else LFit=LP->power[0];
  809.         fprintf(myfile,"%9.4g\t",LFit);
  810.  
  811.         /* quadratic fit */
  812.         LFit=0.0;
  813.         VV=1.0;
  814.         for(j=0;j<3;j++){
  815.             LFit += LP->q[j]*VV;
  816.             VV *= V[i];
  817.         }
  818.         fprintf(myfile,"%9.4g\t",L[i]-LFit);
  819.         /* polynomial fit */
  820.         LFit=0.0;
  821.         VV=1.0;
  822.         for(j=0;j<LP->coefficients;j++){
  823.             LFit += LP->p[j]*VV;
  824.             VV *= V[i];
  825.         }
  826.         fprintf(myfile,"%9.4g\t",L[i]-LFit);
  827.         /* power law fit */
  828.         LFit=LP->power[1]+LP->power[2]*V[i];
  829.         if(LFit>0.0) LFit=LP->power[0]+pow(LFit,LP->power[3]);
  830.         else LFit=LP->power[0];
  831.         fprintf(myfile,"%9.4g\t",L[i]-LFit);
  832.  
  833.         if(i<nLrgb) for(k=0;k<3;k++) {
  834.             VV=LP->r*nrgb[k][i][0]+LP->g*nrgb[k][i][1]+LP->b*nrgb[k][i][2];
  835.             fprintf(myfile,"%f\t",VV);
  836.         }
  837.         if(i<nLrgb) for(k=0;k<3;k++) fprintf(myfile,"%f\t",Lrgb[k][i]);
  838.         if(i<nLrgb) for(k=0;k<3;k++) fprintf(myfile,"%f\t",LToV(LP,Lrgb[k][i]));
  839.         if(i==0) fprintf(myfile,"%6g\t%6g\t%6g\t%6g",LP->r,LP->g,LP->b
  840.                 ,gain[0]+gain[1]+gain[2]);
  841.     }
  842.     fclose(myfile);
  843. }
  844.  
  845. #if 1
  846. void SaveLuminanceRecord(LuminanceRecord *LP)
  847. /* Create C header file "LuminanceRecord?.h", where ? is the screen number. */
  848. {
  849.     FILE *file;
  850.     int i;
  851.     static char string[100],filename[100],dateString[100];
  852.     long seconds;
  853.     
  854.     LP->VMin=0;
  855.     LP->VMax=(1L<<LP->dacSize)-1;    /* maximum value that can be loaded into DAC */
  856.     LP->LMin=VToL(LP,LP->VMin);        /* min luminance */
  857.     LP->LMax=VToL(LP,LP->VMax);        /* max luminance */
  858.     LP->L.exists=0;
  859.     sprintf(filename,"LuminanceRecord%hd.h",LP->screen);
  860.     GetDateTime((void *)&seconds);
  861.     printf("Please enter the following descriptive information for the %s file.\n",filename);
  862.  
  863.     /* time and date */
  864.     IUDateString(seconds,longDate,(unsigned char *)string);
  865.     p2cstr((unsigned char *)string);
  866.     IUTimeString(seconds,FALSE,(unsigned char *)dateString);
  867.     p2cstr((unsigned char *)dateString);
  868.     sprintf(dateString,"%s %s",dateString,string);
  869.     LP->date=dateString;
  870.     printf("LR.date=\"%s\";\n",LP->date);    
  871.  
  872.     printf("Enter LR.id, i.e. monitor model and serial # (%s):",LP->id);
  873.     gets(string);
  874.     if(strlen(string)>0)strcpy(LP->id=malloc(strlen(string)+1),string);
  875.     printf("%s\n",LP->id);
  876.  
  877.     printf("Enter LR.name, i.e. informal monitor name \n(%s):",LP->name);
  878.     gets(string);
  879.     if(strlen(string)>0)strcpy(LP->name=malloc(strlen(string)+1),string);
  880.     printf("%s\n",LP->name);
  881.  
  882.     printf("Enter LR.notes, i.e. details of calibration: who & how\n(%s):\n",LP->notes);
  883.     gets(string);
  884.     if(strlen(string)>0) strcpy(LP->notes=malloc(strlen(string)+1),string);
  885.     printf("%s\n",LP->notes);
  886.  
  887.     printf("Nominally, the Apple High-Res. Monochrome Monitor has 76 pixels/inch.\n");
  888.     printf("Nominally, the Apple High-Res. RGB Monitor has 69 pixels/inch.\n");
  889.     printf("Enter LR.dpi, i.e. pixels per inch (%.1f):",LP->dpi);
  890.     gets(string);
  891.     sscanf(string,"%lf",&LP->dpi);
  892.     printf("%.1f\n",LP->dpi);
  893.  
  894.     printf("Enter LR.units luminance units (%s):",LP->units);
  895.     gets(string);
  896.     if(strlen(string)>0)strcpy(LP->units=malloc(strlen(string)+1),string);
  897.     printf("%s\n",LP->units);
  898.  
  899.     // Write comment at top of file.
  900.     file=fopen(filename,"w");
  901.     fprintf(file,"/* This %s file is a description of a monitor produced by the CalibrateLuminance program. */\n",filename);
  902.     fprintf(file,"/* This file may be #included in any C program. It simply fills in a LuminanceRecord data structure. */\n");
  903.     fprintf(file,"/* Or use ReadLuminanceRecord.c to read this file at runtime. */\n");
  904.     fprintf(file,"/* The theory is described by Pelli & Zhang (1991). The data structure is defined in Luminance.h.  */\n");
  905.     fprintf(file,"/* CalibrateLuminance and Luminance.h are part of the VideoToolbox software. */\n");
  906.     fprintf(file,"/* Pelli & Zhang (1991) Accurate control of contrast on microcomputer displays. */\n");
  907.     fprintf(file,"/* Vision Research, 31:1337-1360. */\n");
  908.     fprintf(file,"/* The VideoToolbox software is available from: */\n");
  909.     fprintf(file,"/* ftp://mirrors.aol.com/pub/info-mac/dev/lib/video-toolbox */\n");
  910.     fprintf(file,"/* Caution: the screen number used here and in GetScreen Device is NOT the same as */\n");
  911.     fprintf(file,"/* displayed by the Monitors cdev in the Control Panel. */\n");
  912.     fclose(file);
  913.     i=WriteLuminanceRecord(filename,LP,0);// assignNoPrintfExit
  914.     // if(i<0)printf("WARNING: WriteLuminanceRecord returned error %d\n",i);
  915. }
  916. #else
  917. {
  918.     FILE *file;
  919.     int i;
  920.     static char string[100],dateString[100];
  921.     long seconds;
  922.     
  923.     LP->VMin=0;
  924.     LP->VMax=(1L<<LP->dacSize)-1;    /* maximum value that can be loaded into DAC */
  925.     LP->LMin=VToL(LP,LP->VMin);        /* min luminance */
  926.     LP->LMax=VToL(LP,LP->VMax);        /* max luminance */
  927.     sprintf(string,"LuminanceRecord%hd.h",LP->screen);
  928.     GetDateTime((void *)&seconds);
  929.     file=fopen(string,"w");
  930.     SetFileInfo(string,'TEXT','KAHL');
  931.     printf("Please enter the following descriptive information for the %s file.\n",string);
  932.  
  933.     fprintf(file,"/* This %s file is a description of a monitor produced by the CalibrateLuminance program. */\n",string);
  934.     fprintf(file,"/* This file may be #included in any C program. It simply fills in a LuminanceRecord data structure. */\n");
  935.     fprintf(file,"/* Or use ReadLuminanceRecord.c to read this file at runtime. */\n");
  936.     fprintf(file,"/* The theory is described by Pelli & Zhang (1991). The data structure is defined in Luminance.h.  */\n");
  937.     fprintf(file,"/* CalibrateLuminance and Luminance.h are part of the VideoToolbox software. */\n");
  938.     fprintf(file,"/* Pelli & Zhang (1991) Accurate control of contrast on microcomputer displays. */\n");
  939.     fprintf(file,"/* Vision Research, 31:1337-1360. */\n");
  940.     fprintf(file,"/* The VideoToolbox software is available from: */\n");
  941.     fprintf(file,"/* ftp://mirrors.aol.com/pub/info-mac/dev/lib/video-toolbox */\n");
  942.     fprintf(file,"/* Caution: the screen number used here and in GetScreen Device is NOT the same as */\n");
  943.     fprintf(file,"/* displayed by the Monitors cdev in the Control Panel. Sorry. The most obvious difference */\n");
  944.     fprintf(file,"/* is that GetScreenDevice always assigns 0 to the main screen, the one with the menu bar. */\n");
  945.     fprintf(file,"LR.screen=%hd;    /* device=GetScreenDevice(LR.screen); */\n"
  946.         ,LP->screen);
  947.     /* time and date */
  948.     IUDateString(seconds,longDate,(unsigned char *)string);
  949.     p2cstr((unsigned char *)string);
  950.     IUTimeString(seconds,FALSE,(unsigned char *)dateString);
  951.     p2cstr((unsigned char *)dateString);
  952.     sprintf(dateString,"%s %s",dateString,string);
  953.     LP->date=dateString;
  954.     fprintf(file,"LR.date=\"%s\";\n",LP->date);
  955.     printf("LR.date=\"%s\";\n",LP->date);    
  956.  
  957.     printf("Enter LR.id, i.e. monitor model and serial # (%s):",LP->id);
  958.     gets(string);
  959.     if(strlen(string)>0)strcpy(LP->id=malloc(strlen(string)+1),string);
  960.     printf("%s\n",LP->id);
  961.     fprintf(file,"LR.id=\"%s\";\n",LP->id);
  962.  
  963.     printf("Enter LR.name, i.e. informal monitor name \n(%s):",LP->name);
  964.     gets(string);
  965.     if(strlen(string)>0)strcpy(LP->name=malloc(strlen(string)+1),string);
  966.     printf("%s\n",LP->name);
  967.     fprintf(file,"LR.name=\"%s\";\n",LP->name);
  968.  
  969.     printf("Enter LR.notes, i.e. details of calibration: who & how\n(%s):\n",LP->notes);
  970.     gets(string);
  971.     if(strlen(string)>0) strcpy(LP->notes=malloc(strlen(string)+1),string);
  972.     printf("%s\n",LP->notes);
  973.     fprintf(file,"LR.notes=\"%s\";\n",LP->notes);
  974.  
  975.     printf("Nominally, the Apple High-Res. Monochrome Monitor has 76 pixels/inch.\n");
  976.     printf("Nominally, the Apple High-Res. RGB Monitor has 69 pixels/inch.\n");
  977.     printf("Enter LR.dpi, i.e. pixels per inch (%.1f):",LP->dpi);
  978.     gets(string);
  979.     sscanf(string,"%lf",&LP->dpi);
  980.     printf("%.1f\n",LP->dpi);
  981.     fprintf(file,"LR.dpi=%.1f;\t/* pixels per inch */\n",LP->dpi);
  982.  
  983.     fprintf(file,"LR.Hz=%.2f;\t/* frames per second */\n",LP->Hz);
  984.  
  985.     printf("Enter LR.units luminance units (%s):",LP->units);
  986.     gets(string);
  987.     if(strlen(string)>0)strcpy(LP->units=malloc(strlen(string)+1),string);
  988.     printf("%s\n",LP->units);
  989.     fprintf(file,"LR.units=\"%s\";\n",LP->units);
  990.  
  991.     fprintf(file,"/* coefficients of polynomial fit */\n");
  992.     fprintf(file,"LR.coefficients=%ld;    /* # of coefficients in polynomial fit */\n"
  993.         ,LP->coefficients);
  994.     fprintf(file,"/* L(V)=p[0]+p[1]*V+p[2]*V*V+ . . . ±polynomialError */\n");
  995.     for(i=0;i<LP->coefficients;i++)    {
  996.         fprintf(file,"LR.p[%d]=%6g;\n",i,LP->p[i]);
  997.     }
  998.     fprintf(file,"LR.polynomialError=%8.4f;    /* RMS error of fit */\n",LP->polynomialError);
  999.     fprintf(file,"/* coefficients of quadratic fit */\n");
  1000.     fprintf(file,"/* L(V)=q[0]+q[1]*V+q[2]*V*V±quadraticError */\n");
  1001.     for(i=0;i<3;i++)    {
  1002.         fprintf(file,"LR.q[%d]=%6g;\n",i,LP->q[i]);
  1003.     }
  1004.     fprintf(file,"LR.quadraticError=%8.4f;    /* RMS error of fit */\n",LP->quadraticError);
  1005.     fprintf(file,"/* coefficients of power law fit */\n");
  1006.     fprintf(file,"/* L(V)=power[0]+Rectify(power[1]+power[2]*V)^power[3]±powerError */\n");
  1007.     fprintf(file,"/* where Rectify(x)=x if x≥0, and Rectify(x)=0 otherwise */\n");
  1008.     fprintf(file,"/* Pelli & Zhang (1991) Eqs.9&10 use symbols v=V/VMax, alpha=power[0], beta=power[1], kappa=power[2]*255, gamma=power[3] */\n");
  1009.     for(i=0;i<4;i++)    {
  1010.         fprintf(file,"LR.power[%d]=%6g;\n",i,LP->power[i]);
  1011.     }
  1012.     fprintf(file,"LR.powerError=%8.4f;    /* RMS error of fit */\n",LP->powerError);
  1013.     fprintf(file,"/* coefficients of power law fit, with fixed exponent */\n");
  1014.     fprintf(file,"/* L(V)=fixedPower[0]+Rectify(fixedPower[1]+fixedPower[2]*V)^fixedPower[3]±fixedPowerError */\n");
  1015.     for(i=0;i<4;i++)    {
  1016.         fprintf(file,"LR.fixedPower[%d]=%6g;\n",i,LP->fixedPower[i]);
  1017.     }
  1018.     fprintf(file,"LR.fixedPowerError=%8.4f;    /* RMS error of fit */\n",LP->fixedPowerError);
  1019.     fprintf(file,"LR.r=%6g;\n",LP->r);
  1020.     fprintf(file,"LR.g=%6g;\n",LP->g);
  1021.     fprintf(file,"LR.b=%6g;\n",LP->b);
  1022.     fprintf(file,"LR.gainAccuracy=%6g;\n",LP->gainAccuracy);
  1023.     fprintf(file,"LR.gm=%6g;    /* The monitor's contrast gain. */\n",LP->gm);
  1024.     fprintf(file,"LR.dacSize=%hd;    /* bits */\n",LP->dacSize);
  1025.     fprintf(file,"LR.VMin=%3d;    /* minimum value that can be loaded into DAC */\n",LP->VMin);
  1026.     fprintf(file,"LR.VMax=%3d;    /* maximum value that can be loaded into DAC */\n",LP->VMax);
  1027.     fprintf(file,"LR.LMin=%8.2f;    /* luminance at VMin */\n",LP->LMin);
  1028.     fprintf(file,"LR.LMax=%8.2f;    /* luminance at VMax */\n",LP->LMax);
  1029.     fprintf(file,"LR.LBackground=%8.3f;    /* background luminance during calibration */\n",LP->LBackground);
  1030.     fprintf(file,"LR.VBackground=%hd;    /* background number used during calibration */\n",LP->VBackground);
  1031.     fprintf(file,"LR.rangeSet=0;    /* indicate that range parameters have yet to be set */\n");
  1032.     fprintf(file,"LR.L.exists=0;    /* indicate that luminance table has yet to be initialized */\n");
  1033.     fclose(file);
  1034. }
  1035. #endif
  1036.  
  1037. FLOAT PowerRMSError(FLOAT pp[])
  1038. /*
  1039. Returns average squared error of power law fit.
  1040. pp[] contains the parameters of a power law fit to the
  1041. luminance measurements L[] made at numbers V[].
  1042. */
  1043. {
  1044.     int i;
  1045.     FLOAT f,e,a;
  1046.     extern double L[];
  1047.     extern int V[];
  1048.     extern int nL;
  1049.     extern int variables;
  1050.     extern FLOAT *p;
  1051.     static FLOAT q[5];
  1052.     
  1053.     for(i=1;i<=variables;i++)q[i]=pp[i];    /* copy the variable parameters */
  1054.     for(i=variables+1;i<=4;i++)q[i]=p[i];    /* copy the fixed parameters */
  1055.     e=0.0;
  1056.     for(i=0;i<nL;i++){
  1057.         f=q[2]+q[3]*V[i];
  1058.         if(f>0.0) f=q[1]+pow(f,q[4]);
  1059.         else f=q[1];
  1060.         a=f-L[i];
  1061.         e += a*a;
  1062.     }
  1063.     e=sqrt(e/nL);
  1064.     return e;
  1065. }
  1066.  
  1067.